!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!

! **************************************************************************************************
!> \brief Routines for optimizing load balance between processes in HFX calculations
!> \par History
!>      04.2008 created [Manuel Guidon]
!>      11.2019 fixed initial value for potential_id (A. Bussy)
!> \author Manuel Guidon
! **************************************************************************************************
MODULE hfx_pair_list_methods
   USE cell_types,                      ONLY: cell_type,&
                                              pbc
   USE gamma,                           ONLY: fgamma => fgamma_0
   USE hfx_types,                       ONLY: &
        hfx_basis_type, hfx_block_range_type, hfx_cell_type, hfx_pgf_list, hfx_pgf_product_list, &
        hfx_potential_type, hfx_screen_coeff_type, pair_list_type, pair_set_list_type
   USE input_constants,                 ONLY: &
        do_potential_TShPSC, do_potential_coulomb, do_potential_gaussian, do_potential_id, &
        do_potential_long, do_potential_mix_cl, do_potential_mix_cl_trunc, do_potential_mix_lg, &
        do_potential_short, do_potential_truncated
   USE kinds,                           ONLY: dp
   USE libint_wrapper,                  ONLY: prim_data_f_size
   USE mathconstants,                   ONLY: pi
   USE mp2_types,                       ONLY: pair_list_type_mp2
   USE particle_types,                  ONLY: particle_type
   USE t_c_g0,                          ONLY: t_c_g0_n
   USE t_sh_p_s_c,                      ONLY: trunc_CS_poly_n20
#include "./base/base_uses.f90"

   IMPLICIT NONE
   PRIVATE

   PUBLIC :: build_pair_list, &
             build_pair_list_mp2, &
             build_pair_list_pgf, &
             build_pgf_product_list, &
             build_atomic_pair_list, &
             pgf_product_list_size

   ! an initial estimate for the size of the product list
   INTEGER, SAVE :: pgf_product_list_size = 128

!***

CONTAINS

! **************************************************************************************************
!> \brief ...
!> \param list1 ...
!> \param list2 ...
!> \param product_list ...
!> \param nproducts ...
!> \param log10_pmax ...
!> \param log10_eps_schwarz ...
!> \param neighbor_cells ...
!> \param cell ...
!> \param potential_parameter ...
!> \param m_max ...
!> \param do_periodic ...
! **************************************************************************************************
   SUBROUTINE build_pgf_product_list(list1, list2, product_list, nproducts, &
                                     log10_pmax, log10_eps_schwarz, neighbor_cells, &
                                     cell, potential_parameter, m_max, do_periodic)

      TYPE(hfx_pgf_list)                                 :: list1, list2
      TYPE(hfx_pgf_product_list), ALLOCATABLE, &
         DIMENSION(:), INTENT(INOUT)                     :: product_list
      INTEGER, INTENT(OUT)                               :: nproducts
      REAL(dp), INTENT(IN)                               :: log10_pmax, log10_eps_schwarz
      TYPE(hfx_cell_type), DIMENSION(:), POINTER         :: neighbor_cells
      TYPE(cell_type), POINTER                           :: cell
      TYPE(hfx_potential_type)                           :: potential_parameter
      INTEGER, INTENT(IN)                                :: m_max
      LOGICAL, INTENT(IN)                                :: do_periodic

      INTEGER                                            :: i, j, k, l, nimages1, nimages2, tmp_i4
      LOGICAL                                            :: use_gamma
      REAL(dp) :: C11(3), den, Eta, EtaInv, factor, Fm(prim_data_f_size), G(3), num, omega2, &
         omega_corr, omega_corr2, P(3), pgf_max_1, pgf_max_2, PQ(3), Q(3), R, R1, R2, ra(3), &
         rb(3), rc(3), rd(3), Rho, RhoInv, rpq2, S1234, S1234a, S1234b, shift(3), ssss, T, &
         temp(3), temp_CC(3), temp_DD(3), tmp, tmp_D(3), W(3), Zeta1, Zeta_C, Zeta_D, ZetapEtaInv
      TYPE(hfx_pgf_product_list), ALLOCATABLE, &
         DIMENSION(:)                                    :: tmp_product_list

      nimages1 = list1%nimages
      nimages2 = list2%nimages
      nproducts = 0
      Zeta1 = list1%zetapzetb
      Eta = list2%zetapzetb
      EtaInv = list2%ZetaInv
      Zeta_C = list2%zeta
      Zeta_D = list2%zetb
      temp_CC = 0.0_dp
      temp_DD = 0.0_dp
      DO i = 1, nimages1
         P = list1%image_list(i)%P
         R1 = list1%image_list(i)%R
         S1234a = list1%image_list(i)%S1234
         pgf_max_1 = list1%image_list(i)%pgf_max
         ra = list1%image_list(i)%ra
         rb = list1%image_list(i)%rb
         DO j = 1, nimages2
            pgf_max_2 = list2%image_list(j)%pgf_max
            IF (pgf_max_1 + pgf_max_2 + log10_pmax < log10_eps_schwarz) CYCLE
            Q = list2%image_list(j)%P
            R2 = list2%image_list(j)%R
            S1234b = list2%image_list(j)%S1234
            rc = list2%image_list(j)%ra
            rd = list2%image_list(j)%rb

            ZetapEtaInv = Zeta1 + Eta
            ZetapEtaInv = 1.0_dp/ZetapEtaInv
            Rho = Zeta1*Eta*ZetapEtaInv
            RhoInv = 1.0_dp/Rho
            S1234 = EXP(S1234a + S1234b)
            IF (do_periodic) THEN
               temp = P - Q
               PQ = pbc(temp, cell)
               shift = -PQ + temp
               temp_CC = rc + shift
               temp_DD = rd + shift
            END IF

            DO k = 1, SIZE(neighbor_cells)
               IF (do_periodic) THEN
                  C11 = temp_CC + neighbor_cells(k)%cell_r(:)
                  tmp_D = temp_DD + neighbor_cells(k)%cell_r(:)
               ELSE
                  C11 = rc
                  tmp_D = rd
               END IF
               Q = (Zeta_C*C11 + Zeta_D*tmp_D)*EtaInv
               rpq2 = (P(1) - Q(1))**2 + (P(2) - Q(2))**2 + (P(3) - Q(3))**2
               IF (potential_parameter%potential_type == do_potential_truncated .OR. &
                   potential_parameter%potential_type == do_potential_short .OR. &
                   potential_parameter%potential_type == do_potential_mix_cl_trunc) THEN
                  IF (rpq2 > (R1 + R2 + potential_parameter%cutoff_radius)**2) CYCLE
               END IF
               IF (potential_parameter%potential_type == do_potential_TShPSC) THEN
                  IF (rpq2 > (R1 + R2 + potential_parameter%cutoff_radius*2.0_dp)**2) CYCLE
               END IF
               nproducts = nproducts + 1

               ! allocate size as needed,
               ! updating the global size estimate to make this a rare event in longer simulations
               IF (nproducts > SIZE(product_list)) THEN
!$OMP              ATOMIC READ
                  tmp_i4 = pgf_product_list_size
                  tmp_i4 = MAX(pgf_product_list_size, (3*nproducts + 1)/2)
!$OMP              ATOMIC WRITE
                  pgf_product_list_size = tmp_i4
                  ALLOCATE (tmp_product_list(SIZE(product_list)))
                  tmp_product_list(:) = product_list
                  DEALLOCATE (product_list)
                  ALLOCATE (product_list(tmp_i4))
                  product_list(1:SIZE(tmp_product_list)) = tmp_product_list
                  DEALLOCATE (tmp_product_list)
               END IF

               T = Rho*rpq2
               SELECT CASE (potential_parameter%potential_type)
               CASE (do_potential_truncated)
                  R = potential_parameter%cutoff_radius*SQRT(Rho)
                  CALL t_c_g0_n(product_list(nproducts)%Fm(1), use_gamma, R, T, m_max)
                  IF (use_gamma) CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
                  factor = 2.0_dp*Pi*RhoInv
               CASE (do_potential_TShPSC)
                  R = potential_parameter%cutoff_radius*SQRT(Rho)
                  product_list(nproducts)%Fm = 0.0_dp
                  CALL trunc_CS_poly_n20(product_list(nproducts)%Fm(1), R, T, m_max)
                  factor = 2.0_dp*Pi*RhoInv
               CASE (do_potential_coulomb)
                  CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
                  factor = 2.0_dp*Pi*RhoInv
               CASE (do_potential_short)
                  CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
                  omega2 = potential_parameter%omega**2
                  omega_corr2 = omega2/(omega2 + Rho)
                  omega_corr = SQRT(omega_corr2)
                  T = T*omega_corr2
                  CALL fgamma(m_max, T, Fm)
                  tmp = -omega_corr
                  DO l = 1, m_max + 1
                     product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l) + Fm(l)*tmp
                     tmp = tmp*omega_corr2
                  END DO
                  factor = 2.0_dp*Pi*RhoInv
               CASE (do_potential_long)
                  omega2 = potential_parameter%omega**2
                  omega_corr2 = omega2/(omega2 + Rho)
                  omega_corr = SQRT(omega_corr2)
                  T = T*omega_corr2
                  CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
                  tmp = omega_corr
                  DO l = 1, m_max + 1
                     product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l)*tmp
                     tmp = tmp*omega_corr2
                  END DO
                  factor = 2.0_dp*Pi*RhoInv
               CASE (do_potential_mix_cl)
                  CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))
                  omega2 = potential_parameter%omega**2
                  omega_corr2 = omega2/(omega2 + Rho)
                  omega_corr = SQRT(omega_corr2)
                  T = T*omega_corr2
                  CALL fgamma(m_max, T, Fm)
                  tmp = omega_corr
                  DO l = 1, m_max + 1
                     product_list(nproducts)%Fm(l) = &
                        product_list(nproducts)%Fm(l)*potential_parameter%scale_coulomb &
                        + Fm(l)*tmp*potential_parameter%scale_longrange
                     tmp = tmp*omega_corr2
                  END DO
                  factor = 2.0_dp*Pi*RhoInv
               CASE (do_potential_mix_cl_trunc)

                  ! truncated
                  R = potential_parameter%cutoff_radius*SQRT(rho)
                  CALL t_c_g0_n(product_list(nproducts)%Fm(1), use_gamma, R, T, m_max)
                  IF (use_gamma) CALL fgamma(m_max, T, product_list(nproducts)%Fm(1))

                  ! Coulomb
                  CALL fgamma(m_max, T, Fm)

                  DO l = 1, m_max + 1
                     product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l)* &
                                                     (potential_parameter%scale_coulomb + potential_parameter%scale_longrange) - &
                                                     Fm(l)*potential_parameter%scale_longrange
                  END DO

                  ! longrange
                  omega2 = potential_parameter%omega**2
                  omega_corr2 = omega2/(omega2 + Rho)
                  omega_corr = SQRT(omega_corr2)
                  T = T*omega_corr2
                  CALL fgamma(m_max, T, Fm)
                  tmp = omega_corr
                  DO l = 1, m_max + 1
                     product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l) + Fm(l)*tmp*potential_parameter%scale_longrange
                     tmp = tmp*omega_corr2
                  END DO
                  factor = 2.0_dp*Pi*RhoInv

               CASE (do_potential_gaussian)
                  omega2 = potential_parameter%omega**2
                  T = -omega2*T/(Rho + omega2)
                  tmp = 1.0_dp
                  DO l = 1, m_max + 1
                     product_list(nproducts)%Fm(l) = EXP(T)*tmp
                     tmp = tmp*omega2/(Rho + omega2)
                  END DO
                  factor = (Pi/(Rho + omega2))**(1.5_dp)
               CASE (do_potential_mix_lg)
                  omega2 = potential_parameter%omega**2
                  omega_corr2 = omega2/(omega2 + Rho)
                  omega_corr = SQRT(omega_corr2)
                  T = T*omega_corr2
                  CALL fgamma(m_max, T, Fm)
                  tmp = omega_corr*2.0_dp*Pi*RhoInv*potential_parameter%scale_longrange
                  DO l = 1, m_max + 1
                     Fm(l) = Fm(l)*tmp
                     tmp = tmp*omega_corr2
                  END DO
                  T = Rho*rpq2
                  T = -omega2*T/(Rho + omega2)
                  tmp = (Pi/(Rho + omega2))**(1.5_dp)*potential_parameter%scale_gaussian
                  DO l = 1, m_max + 1
                     product_list(nproducts)%Fm(l) = EXP(T)*tmp + Fm(l)
                     tmp = tmp*omega2/(Rho + omega2)
                  END DO
                  factor = 1.0_dp
               CASE (do_potential_id)
                  num = list1%zeta*list1%zetb
                  den = list1%zeta + list1%zetb
                  ssss = -num/den*SUM((ra - rb)**2)

                  num = den*Zeta_C
                  den = den + Zeta_C
                  ssss = ssss - num/den*SUM((P - rc)**2)

                  G(:) = (list1%zeta*ra(:) + list1%zetb*rb(:) + Zeta_C*rc(:))/den
                  num = den*Zeta_D
                  den = den + Zeta_D
                  ssss = ssss - num/den*SUM((G - rd)**2)

                  product_list(nproducts)%Fm(:) = EXP(ssss)
                  factor = 1.0_dp
                  IF (S1234 > EPSILON(0.0_dp)) factor = 1.0_dp/S1234
               END SELECT

               tmp = (Pi*ZetapEtaInv)**3
               factor = factor*S1234*SQRT(tmp)

               DO l = 1, m_max + 1
                  product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l)*factor
               END DO

               W = (Zeta1*P + Eta*Q)*ZetapEtaInv
               product_list(nproducts)%ra = ra
               product_list(nproducts)%rb = rb
               product_list(nproducts)%rc = C11
               product_list(nproducts)%rd = tmp_D
               product_list(nproducts)%ZetapEtaInv = ZetapEtaInv
               product_list(nproducts)%Rho = Rho
               product_list(nproducts)%RhoInv = RhoInv
               product_list(nproducts)%P = P
               product_list(nproducts)%Q = Q
               product_list(nproducts)%W = W
               product_list(nproducts)%AB = ra - rb
               product_list(nproducts)%CD = C11 - tmp_D
            END DO
         END DO
      END DO

   END SUBROUTINE build_pgf_product_list

! **************************************************************************************************
!> \brief ...
!> \param npgfa ...
!> \param npgfb ...
!> \param list ...
!> \param zeta ...
!> \param zetb ...
!> \param screen1 ...
!> \param screen2 ...
!> \param pgf ...
!> \param R_pgf ...
!> \param log10_pmax ...
!> \param log10_eps_schwarz ...
!> \param ra ...
!> \param rb ...
!> \param nelements ...
!> \param neighbor_cells ...
!> \param nimages ...
!> \param do_periodic ...
! **************************************************************************************************
   SUBROUTINE build_pair_list_pgf(npgfa, npgfb, list, zeta, zetb, screen1, screen2, pgf, R_pgf, &
                                  log10_pmax, log10_eps_schwarz, ra, rb, nelements, &
                                  neighbor_cells, nimages, do_periodic)
      INTEGER, INTENT(IN)                                :: npgfa, npgfb
      TYPE(hfx_pgf_list), DIMENSION(npgfa*npgfb)         :: list
      REAL(dp), DIMENSION(1:npgfa), INTENT(IN)           :: zeta
      REAL(dp), DIMENSION(1:npgfb), INTENT(IN)           :: zetb
      REAL(dp), INTENT(IN)                               :: screen1(2), screen2(2)
      TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
         POINTER                                         :: pgf, R_pgf
      REAL(dp), INTENT(IN)                               :: log10_pmax, log10_eps_schwarz, ra(3), &
                                                            rb(3)
      INTEGER, INTENT(OUT)                               :: nelements
      TYPE(hfx_cell_type), DIMENSION(:), POINTER         :: neighbor_cells
      INTEGER                                            :: nimages(npgfa*npgfb)
      LOGICAL, INTENT(IN)                                :: do_periodic

      INTEGER                                            :: element_counter, i, ipgf, j, jpgf
      REAL(dp)                                           :: AB(3), im_B(3), pgf_max, rab2, Zeta1, &
                                                            Zeta_A, Zeta_B, ZetaInv

      nimages = 0
      ! ** inner loop may never be reached
      nelements = npgfa*npgfb
      DO i = 1, SIZE(neighbor_cells)
         IF (do_periodic) THEN
            im_B = rb + neighbor_cells(i)%cell_r(:)
         ELSE
            im_B = rb
         END IF
         AB = ra - im_B
         rab2 = AB(1)**2 + AB(2)**2 + AB(3)**2
         IF (screen1(1)*rab2 + screen1(2) + screen2(2) + log10_pmax < log10_eps_schwarz) CYCLE
         element_counter = 0
         DO ipgf = 1, npgfa
            DO jpgf = 1, npgfb
               element_counter = element_counter + 1
               pgf_max = pgf(jpgf, ipgf)%x(1)*rab2 + pgf(jpgf, ipgf)%x(2)
               IF (pgf_max + screen2(2) + log10_pmax < log10_eps_schwarz) THEN
                  CYCLE
               END IF
               nimages(element_counter) = nimages(element_counter) + 1
               list(element_counter)%image_list(nimages(element_counter))%pgf_max = pgf_max
               list(element_counter)%image_list(nimages(element_counter))%ra = ra
               list(element_counter)%image_list(nimages(element_counter))%rb = im_B
               list(element_counter)%image_list(nimages(element_counter))%rab2 = rab2

               Zeta_A = zeta(ipgf)
               Zeta_B = zetb(jpgf)
               Zeta1 = Zeta_A + Zeta_B
               ZetaInv = 1.0_dp/Zeta1

               IF (nimages(element_counter) == 1) THEN
                  list(element_counter)%ipgf = ipgf
                  list(element_counter)%jpgf = jpgf
                  list(element_counter)%zetaInv = ZetaInv
                  list(element_counter)%zetapzetb = Zeta1
                  list(element_counter)%zeta = Zeta_A
                  list(element_counter)%zetb = Zeta_B
               END IF

               list(element_counter)%image_list(nimages(element_counter))%S1234 = (-Zeta_A*Zeta_B*ZetaInv*rab2)
               list(element_counter)%image_list(nimages(element_counter))%P = (Zeta_A*ra + Zeta_B*im_B)*ZetaInv
               list(element_counter)%image_list(nimages(element_counter))%R = &
                  MAX(0.0_dp, R_pgf(jpgf, ipgf)%x(1)*rab2 + R_pgf(jpgf, ipgf)%x(2))
               list(element_counter)%image_list(nimages(element_counter))%ra = ra
               list(element_counter)%image_list(nimages(element_counter))%rb = im_B
               list(element_counter)%image_list(nimages(element_counter))%rab2 = rab2
               list(element_counter)%image_list(nimages(element_counter))%bcell = neighbor_cells(i)%cell
            END DO
         END DO
         nelements = MAX(nelements, element_counter)
      END DO
      DO j = 1, nelements
         list(j)%nimages = nimages(j)
      END DO
      ! ** Remove unused elements

      element_counter = 0
      DO j = 1, nelements
         IF (list(j)%nimages == 0) CYCLE
         element_counter = element_counter + 1
         list(element_counter)%nimages = list(j)%nimages
         list(element_counter)%zetapzetb = list(j)%zetapzetb
         list(element_counter)%ZetaInv = list(j)%ZetaInv
         list(element_counter)%zeta = list(j)%zeta
         list(element_counter)%zetb = list(j)%zetb
         list(element_counter)%ipgf = list(j)%ipgf
         list(element_counter)%jpgf = list(j)%jpgf
         DO i = 1, list(j)%nimages
            list(element_counter)%image_list(i) = list(j)%image_list(i)
         END DO
      END DO

      nelements = element_counter

   END SUBROUTINE build_pair_list_pgf

! **************************************************************************************************
!> \brief ...
!> \param natom ...
!> \param list ...
!> \param set_list ...
!> \param i_start ...
!> \param i_end ...
!> \param j_start ...
!> \param j_end ...
!> \param kind_of ...
!> \param basis_parameter ...
!> \param particle_set ...
!> \param do_periodic ...
!> \param coeffs_set ...
!> \param coeffs_kind ...
!> \param coeffs_kind_max0 ...
!> \param log10_eps_schwarz ...
!> \param cell ...
!> \param pmax_blocks ...
!> \param atomic_pair_list ...
! **************************************************************************************************
   SUBROUTINE build_pair_list(natom, list, set_list, i_start, i_end, j_start, j_end, kind_of, basis_parameter, particle_set, &
                              do_periodic, coeffs_set, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, &
                              pmax_blocks, atomic_pair_list)

      INTEGER, INTENT(IN)                                :: natom
      TYPE(pair_list_type), INTENT(INOUT)                :: list
      TYPE(pair_set_list_type), DIMENSION(:), &
         INTENT(OUT)                                     :: set_list
      INTEGER, INTENT(IN)                                :: i_start, i_end, j_start, j_end, &
                                                            kind_of(*)
      TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
         POINTER                                         :: basis_parameter
      TYPE(particle_type), DIMENSION(:), INTENT(IN), &
         POINTER                                         :: particle_set
      LOGICAL, INTENT(IN)                                :: do_periodic
      TYPE(hfx_screen_coeff_type), &
         DIMENSION(:, :, :, :), INTENT(IN), POINTER      :: coeffs_set
      TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
         INTENT(IN)                                      :: coeffs_kind
      REAL(KIND=dp), INTENT(IN)                          :: coeffs_kind_max0, log10_eps_schwarz
      TYPE(cell_type), POINTER                           :: cell
      REAL(dp), INTENT(IN)                               :: pmax_blocks
      LOGICAL, DIMENSION(natom, natom), INTENT(IN)       :: atomic_pair_list

      INTEGER                                            :: iatom, ikind, iset, jatom, jkind, jset, &
                                                            n_element, nset_ij, nseta, nsetb
      REAL(KIND=dp)                                      :: rab2
      REAL(KIND=dp), DIMENSION(3)                        :: B11, pbc_B, ra, rb, temp

      n_element = 0
      nset_ij = 0

      DO iatom = i_start, i_end
         DO jatom = j_start, j_end
            IF (atomic_pair_list(jatom, iatom) .EQV. .FALSE.) CYCLE

            ikind = kind_of(iatom)
            nseta = basis_parameter(ikind)%nset
            ra = particle_set(iatom)%r(:)

            IF (jatom < iatom) CYCLE
            jkind = kind_of(jatom)
            nsetb = basis_parameter(jkind)%nset
            rb = particle_set(jatom)%r(:)

            IF (do_periodic) THEN
               temp = rb - ra
               pbc_B = pbc(temp, cell)
               B11 = ra + pbc_B
               rab2 = (ra(1) - B11(1))**2 + (ra(2) - B11(2))**2 + (ra(3) - B11(3))**2
            ELSE
               rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2
               B11 = rb ! ra - rb
            END IF
            IF ((coeffs_kind(jkind, ikind)%x(1)*rab2 + &
                 coeffs_kind(jkind, ikind)%x(2)) + coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE

            n_element = n_element + 1
            list%elements(n_element)%pair = [iatom, jatom]
            list%elements(n_element)%kind_pair = [ikind, jkind]
            list%elements(n_element)%r1 = ra
            list%elements(n_element)%r2 = B11
            list%elements(n_element)%dist2 = rab2
            ! build a list of guaranteed overlapping sets
            list%elements(n_element)%set_bounds(1) = nset_ij + 1
            DO iset = 1, nseta
               DO jset = 1, nsetb
                  IF (coeffs_set(jset, iset, jkind, ikind)%x(1)*rab2 + coeffs_set(jset, iset, jkind, ikind)%x(2) + &
                      coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE
                  nset_ij = nset_ij + 1
                  set_list(nset_ij)%pair = [iset, jset]
               END DO
            END DO
            list%elements(n_element)%set_bounds(2) = nset_ij
         END DO
      END DO

      list%n_element = n_element

   END SUBROUTINE build_pair_list

! **************************************************************************************************
!> \brief ...
!> \param natom ...
!> \param atomic_pair_list ...
!> \param kind_of ...
!> \param basis_parameter ...
!> \param particle_set ...
!> \param do_periodic ...
!> \param coeffs_kind ...
!> \param coeffs_kind_max0 ...
!> \param log10_eps_schwarz ...
!> \param cell ...
!> \param blocks ...
! **************************************************************************************************
   SUBROUTINE build_atomic_pair_list(natom, atomic_pair_list, kind_of, basis_parameter, particle_set, &
                                     do_periodic, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, &
                                     blocks)
      INTEGER, INTENT(IN)                                :: natom
      LOGICAL, DIMENSION(natom, natom)                   :: atomic_pair_list
      INTEGER, INTENT(IN)                                :: kind_of(*)
      TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
         POINTER                                         :: basis_parameter
      TYPE(particle_type), DIMENSION(:), INTENT(IN), &
         POINTER                                         :: particle_set
      LOGICAL, INTENT(IN)                                :: do_periodic
      TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
         INTENT(IN)                                      :: coeffs_kind
      REAL(KIND=dp), INTENT(IN)                          :: coeffs_kind_max0, log10_eps_schwarz
      TYPE(cell_type), POINTER                           :: cell
      TYPE(hfx_block_range_type), DIMENSION(:), &
         INTENT(IN), POINTER                             :: blocks

      INTEGER                                            :: iatom, iatom_end, iatom_start, iblock, &
                                                            ikind, jatom, jatom_end, jatom_start, &
                                                            jblock, jkind, nseta, nsetb
      REAL(KIND=dp)                                      :: rab2
      REAL(KIND=dp), DIMENSION(3)                        :: B11, pbc_B, ra, rb, temp

      atomic_pair_list = .FALSE.

      DO iblock = 1, SIZE(blocks)
         iatom_start = blocks(iblock)%istart
         iatom_end = blocks(iblock)%iend
         DO jblock = 1, SIZE(blocks)
            jatom_start = blocks(jblock)%istart
            jatom_end = blocks(jblock)%iend

            DO iatom = iatom_start, iatom_end
               ikind = kind_of(iatom)
               nseta = basis_parameter(ikind)%nset
               ra = particle_set(iatom)%r(:)
               DO jatom = jatom_start, jatom_end
                  IF (jatom < iatom) CYCLE
                  jkind = kind_of(jatom)
                  nsetb = basis_parameter(jkind)%nset
                  rb = particle_set(jatom)%r(:)

                  IF (do_periodic) THEN
                     temp = rb - ra
                     pbc_B = pbc(temp, cell)
                     B11 = ra + pbc_B
                     rab2 = (ra(1) - B11(1))**2 + (ra(2) - B11(2))**2 + (ra(3) - B11(3))**2
                  ELSE
                     rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2
                     B11 = rb ! ra - rb
                  END IF
                  IF ((coeffs_kind(jkind, ikind)%x(1)*rab2 + &
                       coeffs_kind(jkind, ikind)%x(2)) + coeffs_kind_max0 < log10_eps_schwarz) CYCLE

                  atomic_pair_list(jatom, iatom) = .TRUE.
                  atomic_pair_list(iatom, jatom) = .TRUE.
               END DO
            END DO
         END DO
      END DO

   END SUBROUTINE build_atomic_pair_list

! **************************************************************************************************
!> \brief ...
!> \param natom ...
!> \param list ...
!> \param set_list ...
!> \param i_start ...
!> \param i_end ...
!> \param j_start ...
!> \param j_end ...
!> \param kind_of ...
!> \param basis_parameter ...
!> \param particle_set ...
!> \param do_periodic ...
!> \param coeffs_set ...
!> \param coeffs_kind ...
!> \param coeffs_kind_max0 ...
!> \param log10_eps_schwarz ...
!> \param cell ...
!> \param pmax_blocks ...
!> \param atomic_pair_list ...
!> \param skip_atom_symmetry ...
! **************************************************************************************************
   SUBROUTINE build_pair_list_mp2(natom, list, set_list, i_start, i_end, j_start, j_end, kind_of, basis_parameter, particle_set, &
                                  do_periodic, coeffs_set, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, &
                                  pmax_blocks, atomic_pair_list, skip_atom_symmetry)

      INTEGER, INTENT(IN)                                :: natom
      TYPE(pair_list_type_mp2), INTENT(INOUT)            :: list
      TYPE(pair_set_list_type), DIMENSION(:), &
         INTENT(OUT)                                     :: set_list
      INTEGER, INTENT(IN)                                :: i_start, i_end, j_start, j_end, &
                                                            kind_of(*)
      TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
         POINTER                                         :: basis_parameter
      TYPE(particle_type), DIMENSION(:), INTENT(IN), &
         POINTER                                         :: particle_set
      LOGICAL, INTENT(IN)                                :: do_periodic
      TYPE(hfx_screen_coeff_type), &
         DIMENSION(:, :, :, :), INTENT(IN), POINTER      :: coeffs_set
      TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
         INTENT(IN)                                      :: coeffs_kind
      REAL(KIND=dp), INTENT(IN)                          :: coeffs_kind_max0, log10_eps_schwarz
      TYPE(cell_type), POINTER                           :: cell
      REAL(dp), INTENT(IN)                               :: pmax_blocks
      LOGICAL, DIMENSION(natom, natom), INTENT(IN)       :: atomic_pair_list
      LOGICAL, INTENT(IN), OPTIONAL                      :: skip_atom_symmetry

      INTEGER                                            :: iatom, ikind, iset, jatom, jkind, jset, &
                                                            n_element, nset_ij, nseta, nsetb
      LOGICAL                                            :: my_skip_atom_symmetry
      REAL(KIND=dp)                                      :: rab2
      REAL(KIND=dp), DIMENSION(3)                        :: B11, pbc_B, ra, rb, temp

      n_element = 0
      nset_ij = 0

      my_skip_atom_symmetry = .FALSE.
      IF (PRESENT(skip_atom_symmetry)) my_skip_atom_symmetry = skip_atom_symmetry

      DO iatom = i_start, i_end
         DO jatom = j_start, j_end
            IF (atomic_pair_list(jatom, iatom) .EQV. .FALSE.) CYCLE

            ikind = kind_of(iatom)
            nseta = basis_parameter(ikind)%nset
            ra = particle_set(iatom)%r(:)

            IF (jatom < iatom .AND. (.NOT. my_skip_atom_symmetry)) CYCLE
            jkind = kind_of(jatom)
            nsetb = basis_parameter(jkind)%nset
            rb = particle_set(jatom)%r(:)

            IF (do_periodic) THEN
               temp = rb - ra
               pbc_B = pbc(temp, cell)
               B11 = ra + pbc_B
               rab2 = (ra(1) - B11(1))**2 + (ra(2) - B11(2))**2 + (ra(3) - B11(3))**2
            ELSE
               rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2
               B11 = rb ! ra - rb
            END IF
            IF ((coeffs_kind(jkind, ikind)%x(1)*rab2 + &
                 coeffs_kind(jkind, ikind)%x(2)) + coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE

            n_element = n_element + 1
            list%elements(n_element)%pair = [iatom, jatom]
            list%elements(n_element)%kind_pair = [ikind, jkind]
            list%elements(n_element)%r1 = ra
            list%elements(n_element)%r2 = B11
            list%elements(n_element)%dist2 = rab2
            ! build a list of guaranteed overlapping sets
            list%elements(n_element)%set_bounds(1) = nset_ij + 1
            DO iset = 1, nseta
               DO jset = 1, nsetb
                  IF (coeffs_set(jset, iset, jkind, ikind)%x(1)*rab2 + coeffs_set(jset, iset, jkind, ikind)%x(2) + &
                      coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE
                  nset_ij = nset_ij + 1
                  set_list(nset_ij)%pair = [iset, jset]
               END DO
            END DO
            list%elements(n_element)%set_bounds(2) = nset_ij
         END DO
      END DO

      list%n_element = n_element

   END SUBROUTINE build_pair_list_mp2

END MODULE hfx_pair_list_methods
